library(sp)
library(raster)
library(rgeos)
library(sf)
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(leaflet)
library(rnaturalearth)
# load 10 km grid
pg = readRDS("~/Documents/GitHub/purple-air-infiltration/data/grid.RDS")
pg = gBuffer(pg, byid = T, width = 5000, capStyle = "SQUARE")
pg = st_as_sf(pg) %>%
st_transform(4326) %>%
select(ID)
# load land use class
nlcd = read.csv("~/BurkeLab Dropbox/Projects/smoke_PM_prediction/data/2_from_EE/NLCD_areas_10km_grid.csv") %>%
select(grid_id = ID, groups) %>%
mutate(groups = gsub("\\[|\\]", "", groups), # get rid of outer brackets
groups = strsplit(groups, "\\}, \\{")) %>% # split on commas between brackets
unnest(groups, keep_empty = TRUE) %>% # groups is now a list that we want to unnest (i.e. lengthen)
mutate(groups = gsub("\\{|\\}", "", groups)) %>% # drop the extra brackets left behind
separate(groups, into = c("landcover", "area"), sep = ",") %>% # split in commas to get land cover class and area
mutate(landcover = trimws(gsub("landcover=", "", landcover, fixed = TRUE)), # drop "landcover"
area = trimws(gsub("sum=", "", area, fixed = TRUE)) %>% as.numeric, # drop "sum"
landcover = recode(landcover, # recode the landcover variables to their classes
"1.0" = "water",
"2.0" = "developed",
"3.0" = "barren",
"4.0" = "forest",
"5.0" = "shrubland",
"7.0" = "herbaceous",
"8.0" = "cultivated",
"9.0" = "wetlands")) %>%
pivot_wider(names_from = landcover, values_from = area, # make it wider, one row for each grid cell, filling missings with 0s because that land class wasn't in the grid cell
values_fill = 0) %>%
mutate(total = water + developed + barren + forest + shrubland + herbaceous + cultivated + wetlands) %>% # calculate total area for the grid cell
mutate(across(!total & !grid_id, ~.x/total)) # calculate percentages in each landcover class
# get land use class as sf object
pg_nlcd = full_join(pg, nlcd %>% select(ID = grid_id, nlcd = barren), by = "ID")
# plot and color 10 km grid cells missing land use class
ggplot(pg_nlcd) + geom_sf(aes(fill = is.na(nlcd)), color = NA)

729 (0.73%) of our 10 km grid cells are missing NLCD. They are located throughout the border and coast.
# load elevation
elevation = read.csv("~/BurkeLab Dropbox/Projects/smoke_PM_prediction/data/2_from_EE/elevation_avg_sd_10km_grid.csv")
# get elevation as sf object
pg_elevation = full_join(pg, elevation %>% select(ID, elevation = mean), by = "ID")
# plot and color 10 km grid cells missing elevation
ggplot(pg_elevation) + geom_sf(aes(fill = is.na(elevation)), color = NA)

128 (0.13%) of our 10 km grid cells are missing elevation. They are located at the Canadian border, the Gulf of Mexico, and the Atlantic coast.
# load era5-land
era5l = readRDS("~/BurkeLab Dropbox/Data/ERA5/Land/2m_temperature/USA/10km_grid/UTC-0600/daily_mean_of_1-hourly/grid_2m_temperature_daily_mean_2006_01.rds") %>%
filter(date == min(date))
# get era5-land as sf object
pg_era5l = full_join(pg, era5l %>% select(ID = id_grid, era5l = `2m_temperature`), by = "ID")
# plot and color by 10 km grid cells missing era5-land
ggplot(pg_era5l) + geom_sf(aes(fill = is.na(era5l)), color = NA)

31 (0.03%) of our 10 km grid cells are missing ERA5-Land (after matching to NN non-NA within 1 degree). They are located in the Florida Keys.
# get 10 km grid cells missing any of these inputs
pg_na = reduce(list(st_drop_geometry(pg_nlcd),
st_drop_geometry(pg_elevation),
st_drop_geometry(pg_era5l)),
full_join,
by = "ID")
pg_na = full_join(pg, pg_na)
pg_na = pg_na %>% filter(is.na(nlcd) | is.na(elevation) | is.na(era5l))
nrow(pg_na)
[1] 788
In total, 788 of our 10 km grid cells are missing NLCD, elevation, or ERA5-Land.
# load EPA stations
epa = read_sf("~/BurkeLab Dropbox/Data/PM25/epa_station_locations/")
# check if any of those NA grid cells have epa stations in them
pg_epa = epa %>% st_filter(pg_na)
nrow(pg_epa) == 0
[1] TRUE
There are no EPA stations located in the 10 km grid cells missing NLCD, elevation, or ERA5-Land.
ggplot() +
geom_sf(data = pg_na, color = "blue") +
geom_sf(data = epa, color = "green", size = 0.3) +
coord_sf(xlim = c(-125, -63), ylim = c(23, 50))

- Blue = 10 km grid cells missing NLCD, elevation, or ERA5-Land
- Green = EPA stations
# load US landmass
us_land = ne_countries(scale = 10,
country = "United States of America",
returnclass = "sf")
# check if any of those NA grid cells overlap US landmass
pg_na_us_land = pg_na %>% st_filter(us_land)
nrow(pg_na_us_land)
[1] 56
56 (7.11%) of the 10 km grid cells missing NLCD, elevation, or ERA5-Land overlap with US landmass.
ggplot() +
geom_sf(data = us_land, fill = "green") +
geom_sf(data = pg_na, color = "blue", fill = NA) +
coord_sf(xlim = c(-125, -63), ylim = c(23, 50))

- Blue = 10 km grid cells missing NLCD, elevation, or ERA5-Land
- Green = US landmass
ggplot() +
geom_sf(data = us_land, fill = "green") +
geom_sf(data = pg_na_us_land, color = "blue", fill = NA) +
coord_sf(xlim = c(-125, -63), ylim = c(23, 50))

Showing only the grid cells that overlap US landmass.
# load population grid
population = readRDS("~/BurkeLab Dropbox/Data/10km_grid_data/grid_population.RDS")
# get population as sf object
pg_population = full_join(pg, population %>% select(ID = id, population = pop), by = "ID")
# check if any of those NA grid cells have population in them
# caveat: population data is from Gridded Population of the World, so may
# include ex-US population too
pg_na_population = pg_population %>% right_join(st_drop_geometry(pg_na), by = "ID")
pg_na_population %>%
st_drop_geometry() %>%
count(population > 0)
282 (35.79%) of the 10 km grid cells missing NLCD, elevation, or ERA5-Land have population in them. (Caveat: this is based on global population data, so if you only want to count US population, this may not be exactly helpful.)
ggplot() +
geom_sf(data = pg_na_population %>% filter(population > 0),
mapping = aes(fill = population),
color = NA)

They are located at various spots along the border.
# interactive look at the most populated places
populated_places = pg_na_population %>%
filter(population > quantile(population, 0.99)) %>%
arrange(desc(population)) %>%
st_centroid() %>%
st_coordinates() %>%
as.data.frame()
leaflet() %>%
addTiles() %>%
addMarkers(lng = populated_places$X, lat = populated_places$Y)
These are the 10 km grid cells missing NLCD, elevation, or ERA5-Land that have the most population.
LS0tCnRpdGxlOiAiQ2hlY2sgR3JpZCBDZWxscyB3aXRoIE1pc3NpbmcgSW5wdXRzIChMYW5kIFVzZSBDbGFzcywgRWxldmF0aW9uLCBhbmQgRVJBNS1MYW5kKSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGUgPSBGfQprbml0cjo6b3B0c19jaHVuayRzZXQobWVzc2FnZSA9IEYpCmtuaXRyOjpvcHRzX2NodW5rJHNldCh3YXJuaW5nID0gRikKYGBgCgpgYGB7ciBsb2FkX3BhY2thZ2VzfQpsaWJyYXJ5KHNwKQpsaWJyYXJ5KHJhc3RlcikKbGlicmFyeShyZ2VvcykKbGlicmFyeShzZikKbGlicmFyeShkcGx5cikKbGlicmFyeSh0aWR5cikKbGlicmFyeShwdXJycikKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGxlYWZsZXQpCmxpYnJhcnkocm5hdHVyYWxlYXJ0aCkKYGBgCgpgYGB7cn0KIyBsb2FkIDEwIGttIGdyaWQKcGcgPSByZWFkUkRTKCJ+L0RvY3VtZW50cy9HaXRIdWIvcHVycGxlLWFpci1pbmZpbHRyYXRpb24vZGF0YS9ncmlkLlJEUyIpCnBnID0gZ0J1ZmZlcihwZywgYnlpZCA9IFQsIHdpZHRoID0gNTAwMCwgY2FwU3R5bGUgPSAiU1FVQVJFIikKcGcgPSBzdF9hc19zZihwZykgJT4lIAogIHN0X3RyYW5zZm9ybSg0MzI2KSAlPiUgCiAgc2VsZWN0KElEKQpgYGAKCmBgYHtyfQojIGxvYWQgbGFuZCB1c2UgY2xhc3MKbmxjZCA9IHJlYWQuY3N2KCJ+L0J1cmtlTGFiIERyb3Bib3gvUHJvamVjdHMvc21va2VfUE1fcHJlZGljdGlvbi9kYXRhLzJfZnJvbV9FRS9OTENEX2FyZWFzXzEwa21fZ3JpZC5jc3YiKSAlPiUgCiAgc2VsZWN0KGdyaWRfaWQgPSBJRCwgZ3JvdXBzKSAlPiUgCiAgbXV0YXRlKGdyb3VwcyA9IGdzdWIoIlxcW3xcXF0iLCAiIiwgZ3JvdXBzKSwgIyBnZXQgcmlkIG9mIG91dGVyIGJyYWNrZXRzCiAgICAgICAgIGdyb3VwcyA9IHN0cnNwbGl0KGdyb3VwcywgIlxcfSwgXFx7IikpICU+JSAjIHNwbGl0IG9uIGNvbW1hcyBiZXR3ZWVuIGJyYWNrZXRzCiAgdW5uZXN0KGdyb3Vwcywga2VlcF9lbXB0eSA9IFRSVUUpICU+JSAjIGdyb3VwcyBpcyBub3cgYSBsaXN0IHRoYXQgd2Ugd2FudCB0byB1bm5lc3QgKGkuZS4gbGVuZ3RoZW4pCiAgbXV0YXRlKGdyb3VwcyA9IGdzdWIoIlxce3xcXH0iLCAiIiwgZ3JvdXBzKSkgJT4lICAjIGRyb3AgdGhlIGV4dHJhIGJyYWNrZXRzIGxlZnQgYmVoaW5kCiAgc2VwYXJhdGUoZ3JvdXBzLCBpbnRvID0gYygibGFuZGNvdmVyIiwgImFyZWEiKSwgc2VwID0gIiwiKSAlPiUgIyBzcGxpdCBpbiBjb21tYXMgdG8gZ2V0IGxhbmQgY292ZXIgY2xhc3MgYW5kIGFyZWEKICBtdXRhdGUobGFuZGNvdmVyID0gdHJpbXdzKGdzdWIoImxhbmRjb3Zlcj0iLCAiIiwgbGFuZGNvdmVyLCBmaXhlZCA9IFRSVUUpKSwgIyBkcm9wICJsYW5kY292ZXIiCiAgICAgICAgIGFyZWEgPSB0cmltd3MoZ3N1Yigic3VtPSIsICIiLCBhcmVhLCBmaXhlZCA9IFRSVUUpKSAlPiUgYXMubnVtZXJpYywgIyBkcm9wICJzdW0iCiAgICAgICAgIGxhbmRjb3ZlciA9IHJlY29kZShsYW5kY292ZXIsICMgcmVjb2RlIHRoZSBsYW5kY292ZXIgdmFyaWFibGVzIHRvIHRoZWlyIGNsYXNzZXMKICAgICAgICAgICAgICAgICAgICAgICAgICAgICIxLjAiID0gIndhdGVyIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICIyLjAiID0gImRldmVsb3BlZCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAiMy4wIiA9ICJiYXJyZW4iLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIjQuMCIgPSAiZm9yZXN0IiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICI1LjAiID0gInNocnVibGFuZCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAiNy4wIiA9ICJoZXJiYWNlb3VzIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICI4LjAiID0gImN1bHRpdmF0ZWQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIjkuMCIgPSAid2V0bGFuZHMiKSkgJT4lCiAgcGl2b3Rfd2lkZXIobmFtZXNfZnJvbSA9IGxhbmRjb3ZlciwgdmFsdWVzX2Zyb20gPSBhcmVhLCAjIG1ha2UgaXQgd2lkZXIsIG9uZSByb3cgZm9yIGVhY2ggZ3JpZCBjZWxsLCBmaWxsaW5nIG1pc3NpbmdzIHdpdGggMHMgYmVjYXVzZSB0aGF0IGxhbmQgY2xhc3Mgd2Fzbid0IGluIHRoZSBncmlkIGNlbGwKICAgICAgICAgICAgICB2YWx1ZXNfZmlsbCA9IDApICU+JQogIG11dGF0ZSh0b3RhbCA9IHdhdGVyICsgZGV2ZWxvcGVkICsgYmFycmVuICsgZm9yZXN0ICsgc2hydWJsYW5kICsgaGVyYmFjZW91cyArIGN1bHRpdmF0ZWQgKyB3ZXRsYW5kcykgJT4lICMgY2FsY3VsYXRlIHRvdGFsIGFyZWEgZm9yIHRoZSBncmlkIGNlbGwKICBtdXRhdGUoYWNyb3NzKCF0b3RhbCAmICFncmlkX2lkLCB+LngvdG90YWwpKSAgIyBjYWxjdWxhdGUgcGVyY2VudGFnZXMgaW4gZWFjaCBsYW5kY292ZXIgY2xhc3MKCiMgZ2V0IGxhbmQgdXNlIGNsYXNzIGFzIHNmIG9iamVjdApwZ19ubGNkID0gZnVsbF9qb2luKHBnLCBubGNkICU+JSBzZWxlY3QoSUQgPSBncmlkX2lkLCBubGNkID0gYmFycmVuKSwgYnkgPSAiSUQiKQpgYGAKCmBgYHtyfQojIHBsb3QgYW5kIGNvbG9yIDEwIGttIGdyaWQgY2VsbHMgbWlzc2luZyBsYW5kIHVzZSBjbGFzcwpnZ3Bsb3QocGdfbmxjZCkgKyBnZW9tX3NmKGFlcyhmaWxsID0gaXMubmEobmxjZCkpLCBjb2xvciA9IE5BKQpgYGAKCmByIHBnX25sY2QgJT4lIGZpbHRlcihpcy5uYShubGNkKSkgJT4lIG5yb3coKWAgKGByIHJvdW5kKChwZ19ubGNkICU+JSBmaWx0ZXIoaXMubmEobmxjZCkpICU+JSBucm93KCkpLyhwZyAlPiUgbnJvdygpKSoxMDAsIDIpYCUpIG9mIG91ciAxMCBrbSBncmlkIGNlbGxzIGFyZSBtaXNzaW5nIE5MQ0QuIFRoZXkgYXJlIGxvY2F0ZWQgdGhyb3VnaG91dCB0aGUgYm9yZGVyIGFuZCBjb2FzdC4KCmBgYHtyfQojIGxvYWQgZWxldmF0aW9uCmVsZXZhdGlvbiA9IHJlYWQuY3N2KCJ+L0J1cmtlTGFiIERyb3Bib3gvUHJvamVjdHMvc21va2VfUE1fcHJlZGljdGlvbi9kYXRhLzJfZnJvbV9FRS9lbGV2YXRpb25fYXZnX3NkXzEwa21fZ3JpZC5jc3YiKQoKIyBnZXQgZWxldmF0aW9uIGFzIHNmIG9iamVjdApwZ19lbGV2YXRpb24gPSBmdWxsX2pvaW4ocGcsIGVsZXZhdGlvbiAlPiUgc2VsZWN0KElELCBlbGV2YXRpb24gPSBtZWFuKSwgYnkgPSAiSUQiKQpgYGAKCmBgYHtyfQojIHBsb3QgYW5kIGNvbG9yIDEwIGttIGdyaWQgY2VsbHMgbWlzc2luZyBlbGV2YXRpb24KZ2dwbG90KHBnX2VsZXZhdGlvbikgKyBnZW9tX3NmKGFlcyhmaWxsID0gaXMubmEoZWxldmF0aW9uKSksIGNvbG9yID0gTkEpCmBgYAoKYHIgcGdfZWxldmF0aW9uICU+JSBmaWx0ZXIoaXMubmEoZWxldmF0aW9uKSkgJT4lIG5yb3coKWAgKGByIHJvdW5kKChwZ19lbGV2YXRpb24gJT4lIGZpbHRlcihpcy5uYShlbGV2YXRpb24pKSAlPiUgbnJvdygpKS8ocGcgJT4lIG5yb3coKSkqMTAwLCAyKWAlKSBvZiBvdXIgMTAga20gZ3JpZCBjZWxscyBhcmUgbWlzc2luZyBlbGV2YXRpb24uIFRoZXkgYXJlIGxvY2F0ZWQgYXQgdGhlIENhbmFkaWFuIGJvcmRlciwgdGhlIEd1bGYgb2YgTWV4aWNvLCBhbmQgdGhlIEF0bGFudGljIGNvYXN0LgoKYGBge3J9CiMgbG9hZCBlcmE1LWxhbmQKZXJhNWwgPSByZWFkUkRTKCJ+L0J1cmtlTGFiIERyb3Bib3gvRGF0YS9FUkE1L0xhbmQvMm1fdGVtcGVyYXR1cmUvVVNBLzEwa21fZ3JpZC9VVEMtMDYwMC9kYWlseV9tZWFuX29mXzEtaG91cmx5L2dyaWRfMm1fdGVtcGVyYXR1cmVfZGFpbHlfbWVhbl8yMDA2XzAxLnJkcyIpICU+JSAKICBmaWx0ZXIoZGF0ZSA9PSBtaW4oZGF0ZSkpCgojIGdldCBlcmE1LWxhbmQgYXMgc2Ygb2JqZWN0CnBnX2VyYTVsID0gZnVsbF9qb2luKHBnLCBlcmE1bCAlPiUgc2VsZWN0KElEID0gaWRfZ3JpZCwgZXJhNWwgPSBgMm1fdGVtcGVyYXR1cmVgKSwgYnkgPSAiSUQiKQpgYGAKCmBgYHtyfQojIHBsb3QgYW5kIGNvbG9yIGJ5IDEwIGttIGdyaWQgY2VsbHMgbWlzc2luZyBlcmE1LWxhbmQKZ2dwbG90KHBnX2VyYTVsKSArIGdlb21fc2YoYWVzKGZpbGwgPSBpcy5uYShlcmE1bCkpLCBjb2xvciA9IE5BKQpgYGAKCmByIHBnX2VyYTVsICU+JSBmaWx0ZXIoaXMubmEoZXJhNWwpKSAlPiUgbnJvdygpYCAoYHIgcm91bmQoKHBnX2VyYTVsICU+JSBmaWx0ZXIoaXMubmEoZXJhNWwpKSAlPiUgbnJvdygpKS8ocGcgJT4lIG5yb3coKSkqMTAwLCAyKWAlKSBvZiBvdXIgMTAga20gZ3JpZCBjZWxscyBhcmUgbWlzc2luZyBFUkE1LUxhbmQgKGFmdGVyIG1hdGNoaW5nIHRvIE5OIG5vbi1OQSB3aXRoaW4gMSBkZWdyZWUpLiBUaGV5IGFyZSBsb2NhdGVkIGluIHRoZSBGbG9yaWRhIEtleXMuCgpgYGB7cn0KIyBnZXQgMTAga20gZ3JpZCBjZWxscyBtaXNzaW5nIGFueSBvZiB0aGVzZSBpbnB1dHMKcGdfbmEgPSByZWR1Y2UobGlzdChzdF9kcm9wX2dlb21ldHJ5KHBnX25sY2QpLCAKICAgICAgICAgICAgICAgICAgICBzdF9kcm9wX2dlb21ldHJ5KHBnX2VsZXZhdGlvbiksIAogICAgICAgICAgICAgICAgICAgIHN0X2Ryb3BfZ2VvbWV0cnkocGdfZXJhNWwpKSwKICAgICAgICAgICAgICAgZnVsbF9qb2luLCAKICAgICAgICAgICAgICAgYnkgPSAiSUQiKQpwZ19uYSA9IGZ1bGxfam9pbihwZywgcGdfbmEpCnBnX25hID0gcGdfbmEgJT4lIGZpbHRlcihpcy5uYShubGNkKSB8IGlzLm5hKGVsZXZhdGlvbikgfCBpcy5uYShlcmE1bCkpCm5yb3cocGdfbmEpCmBgYAoKSW4gdG90YWwsIGByIHBnX25hICU+JSBucm93KClgIG9mIG91ciAxMCBrbSBncmlkIGNlbGxzIGFyZSBtaXNzaW5nIE5MQ0QsIGVsZXZhdGlvbiwgb3IgRVJBNS1MYW5kLgoKYGBge3J9CiMgbG9hZCBFUEEgc3RhdGlvbnMKZXBhID0gcmVhZF9zZigifi9CdXJrZUxhYiBEcm9wYm94L0RhdGEvUE0yNS9lcGFfc3RhdGlvbl9sb2NhdGlvbnMvIikKCiMgY2hlY2sgaWYgYW55IG9mIHRob3NlIE5BIGdyaWQgY2VsbHMgaGF2ZSBlcGEgc3RhdGlvbnMgaW4gdGhlbQpwZ19lcGEgPSBlcGEgJT4lIHN0X2ZpbHRlcihwZ19uYSkKbnJvdyhwZ19lcGEpID09IDAKYGBgCgpUaGVyZSBhcmUgYHIgaWZlbHNlKChwZ19lcGEgJT4lIG5yb3coKSkgPT0gMCwgIm5vIiwgcGdfZXBhICU+JSBucm93KCkpYCBFUEEgc3RhdGlvbnMgbG9jYXRlZCBpbiB0aGUgMTAga20gZ3JpZCBjZWxscyBtaXNzaW5nIE5MQ0QsIGVsZXZhdGlvbiwgb3IgRVJBNS1MYW5kLgoKYGBge3J9CiMgcGxvdCBOQSBncmlkIGNlbGxzIGFuZCBFUEEgc3RhdGlvbnMKZ2dwbG90KCkgKyAKICBnZW9tX3NmKGRhdGEgPSBwZ19uYSwgY29sb3IgPSAiYmx1ZSIpICsgCiAgZ2VvbV9zZihkYXRhID0gZXBhLCBjb2xvciA9ICJncmVlbiIsIHNpemUgPSAwLjMpICsgCiAgY29vcmRfc2YoeGxpbSA9IGMoLTEyNSwgLTYzKSwgeWxpbSA9IGMoMjMsIDUwKSkKYGBgCgoqICBCbHVlID0gMTAga20gZ3JpZCBjZWxscyBtaXNzaW5nIE5MQ0QsIGVsZXZhdGlvbiwgb3IgRVJBNS1MYW5kCiogIEdyZWVuID0gRVBBIHN0YXRpb25zCgpgYGB7cn0KIyBsb2FkIFVTIGxhbmRtYXNzCnVzX2xhbmQgPSBuZV9jb3VudHJpZXMoc2NhbGUgPSAxMCwKICAgICAgICAgICAgICAgICAgICAgICBjb3VudHJ5ID0gIlVuaXRlZCBTdGF0ZXMgb2YgQW1lcmljYSIsCiAgICAgICAgICAgICAgICAgICAgICAgcmV0dXJuY2xhc3MgPSAic2YiKQoKIyBjaGVjayBpZiBhbnkgb2YgdGhvc2UgTkEgZ3JpZCBjZWxscyBvdmVybGFwIFVTIGxhbmRtYXNzCnBnX25hX3VzX2xhbmQgPSBwZ19uYSAlPiUgc3RfZmlsdGVyKHVzX2xhbmQpCm5yb3cocGdfbmFfdXNfbGFuZCkKYGBgCgpgciBwZ19uYV91c19sYW5kICU+JSBucm93KClgIChgciByb3VuZCgocGdfbmFfdXNfbGFuZCAlPiUgbnJvdygpKS8ocGdfbmEgJT4lIG5yb3coKSkqMTAwLCAyKWAlKSBvZiB0aGUgMTAga20gZ3JpZCBjZWxscyBtaXNzaW5nIE5MQ0QsIGVsZXZhdGlvbiwgb3IgRVJBNS1MYW5kIG92ZXJsYXAgd2l0aCBVUyBsYW5kbWFzcy4KCmBgYHtyfQojIHBsb3QgTkEgZ3JpZCBjZWxscyBhbmQgVVMgbGFuZG1hc3MKZ2dwbG90KCkgKyAKICBnZW9tX3NmKGRhdGEgPSB1c19sYW5kLCBmaWxsID0gImdyZWVuIikgKyAKICBnZW9tX3NmKGRhdGEgPSBwZ19uYSwgY29sb3IgPSAiYmx1ZSIsIGZpbGwgPSBOQSkgKyAKICBjb29yZF9zZih4bGltID0gYygtMTI1LCAtNjMpLCB5bGltID0gYygyMywgNTApKQpgYGAKCiogIEJsdWUgPSAxMCBrbSBncmlkIGNlbGxzIG1pc3NpbmcgTkxDRCwgZWxldmF0aW9uLCBvciBFUkE1LUxhbmQKKiAgR3JlZW4gPSBVUyBsYW5kbWFzcwoKYGBge3J9CiMgcGxvdCBvbmx5IHdoZXJlIHRoZXJlIGlzIG92ZXJsYXAKZ2dwbG90KCkgKyAKICBnZW9tX3NmKGRhdGEgPSB1c19sYW5kLCBmaWxsID0gImdyZWVuIikgKyAKICBnZW9tX3NmKGRhdGEgPSBwZ19uYV91c19sYW5kLCBjb2xvciA9ICJibHVlIiwgZmlsbCA9IE5BKSArIAogIGNvb3JkX3NmKHhsaW0gPSBjKC0xMjUsIC02MyksIHlsaW0gPSBjKDIzLCA1MCkpCmBgYAoKU2hvd2luZyBvbmx5IHRoZSBncmlkIGNlbGxzIHRoYXQgb3ZlcmxhcCBVUyBsYW5kbWFzcy4KCmBgYHtyfQojIGxvYWQgcG9wdWxhdGlvbiBncmlkCnBvcHVsYXRpb24gPSByZWFkUkRTKCJ+L0J1cmtlTGFiIERyb3Bib3gvRGF0YS8xMGttX2dyaWRfZGF0YS9ncmlkX3BvcHVsYXRpb24uUkRTIikKCiMgZ2V0IHBvcHVsYXRpb24gYXMgc2Ygb2JqZWN0CnBnX3BvcHVsYXRpb24gPSBmdWxsX2pvaW4ocGcsIHBvcHVsYXRpb24gJT4lIHNlbGVjdChJRCA9IGlkLCBwb3B1bGF0aW9uID0gcG9wKSwgYnkgPSAiSUQiKQoKIyBjaGVjayBpZiBhbnkgb2YgdGhvc2UgTkEgZ3JpZCBjZWxscyBoYXZlIHBvcHVsYXRpb24gaW4gdGhlbQojIGNhdmVhdDogcG9wdWxhdGlvbiBkYXRhIGlzIGZyb20gR3JpZGRlZCBQb3B1bGF0aW9uIG9mIHRoZSBXb3JsZCwgc28gbWF5IAojIGluY2x1ZGUgZXgtVVMgcG9wdWxhdGlvbiB0b28KcGdfbmFfcG9wdWxhdGlvbiA9IHBnX3BvcHVsYXRpb24gJT4lIHJpZ2h0X2pvaW4oc3RfZHJvcF9nZW9tZXRyeShwZ19uYSksIGJ5ID0gIklEIikKcGdfbmFfcG9wdWxhdGlvbiAlPiUgCiAgc3RfZHJvcF9nZW9tZXRyeSgpICU+JSAKICBjb3VudChwb3B1bGF0aW9uID4gMCkKYGBgCgpgciBwZ19uYV9wb3B1bGF0aW9uICU+JSBmaWx0ZXIocG9wdWxhdGlvbiA+IDApICU+JSBucm93KClgIChgciByb3VuZCgocGdfbmFfcG9wdWxhdGlvbiAlPiUgZmlsdGVyKHBvcHVsYXRpb24gPiAwKSAlPiUgbnJvdygpKS8ocGdfbmFfcG9wdWxhdGlvbiAlPiUgbnJvdygpKSoxMDAsIDIpYCUpIG9mIHRoZSAxMCBrbSBncmlkIGNlbGxzIG1pc3NpbmcgTkxDRCwgZWxldmF0aW9uLCBvciBFUkE1LUxhbmQgaGF2ZSBwb3B1bGF0aW9uIGluIHRoZW0uIChDYXZlYXQ6IHRoaXMgaXMgYmFzZWQgb24gZ2xvYmFsIHBvcHVsYXRpb24gZGF0YSwgc28gaWYgeW91IG9ubHkgd2FudCB0byBjb3VudCBVUyBwb3B1bGF0aW9uLCB0aGlzIG1heSBub3QgYmUgZXhhY3RseSBoZWxwZnVsLikKCmBgYHtyfQojIHBsb3Qgb25seSB3aGVyZSB0aGVyZSBpcyBwb3B1bGF0aW9uCmdncGxvdCgpICsgCiAgZ2VvbV9zZihkYXRhID0gcGdfbmFfcG9wdWxhdGlvbiAlPiUgZmlsdGVyKHBvcHVsYXRpb24gPiAwKSwgCiAgICAgICAgICBtYXBwaW5nID0gYWVzKGZpbGwgPSBwb3B1bGF0aW9uKSwKICAgICAgICAgIGNvbG9yID0gTkEpCmBgYAoKVGhleSBhcmUgbG9jYXRlZCBhdCB2YXJpb3VzIHNwb3RzIGFsb25nIHRoZSBib3JkZXIuCgpgYGB7cn0KIyBpbnRlcmFjdGl2ZSBsb29rIGF0IHRoZSBtb3N0IHBvcHVsYXRlZCBwbGFjZXMKcG9wdWxhdGVkX3BsYWNlcyA9IHBnX25hX3BvcHVsYXRpb24gJT4lIAogIGZpbHRlcihwb3B1bGF0aW9uID4gcXVhbnRpbGUocG9wdWxhdGlvbiwgMC45OSkpICU+JSAKICBhcnJhbmdlKGRlc2MocG9wdWxhdGlvbikpICU+JSAKICBzdF9jZW50cm9pZCgpICU+JSAKICBzdF9jb29yZGluYXRlcygpICU+JSAKICBhcy5kYXRhLmZyYW1lKCkKbGVhZmxldCgpICU+JQogIGFkZFRpbGVzKCkgJT4lIAogIGFkZE1hcmtlcnMobG5nID0gcG9wdWxhdGVkX3BsYWNlcyRYLCBsYXQgPSBwb3B1bGF0ZWRfcGxhY2VzJFkpCmBgYAoKVGhlc2UgYXJlIHRoZSAxMCBrbSBncmlkIGNlbGxzIG1pc3NpbmcgTkxDRCwgZWxldmF0aW9uLCBvciBFUkE1LUxhbmQgdGhhdCBoYXZlIHRoZSBtb3N0IHBvcHVsYXRpb24uCg==